/*
    Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// FITS I/O routines for the grid classes.

// $Id$

#include "grid.h"
#include "CCfits/CCfits"
#include "blitz-fits.h"


/** Saves grid structure to a FITS HDU.  The structure is saved as a
    binary table with a Boolean column indicating whether or not a
    cell is subdivided.  The grid cells are traversed in a depth-first
    fashion.  The quantities of the base grid are saved as keywords.
    See the makegrid-readme for more details. */ 
template<typename cell_data_type/*, typename grid_type*/> 
void mcrx::adaptive_grid<cell_data_type/*,grid_type*/>:: 
save_structure (CCfits::FITS& file, const std::string&hdu,
		const std::string& length_unit) const
{
  using namespace CCfits;
  std::cout << "saving grid structure in HDU " << hdu << std::endl; 
  const string length_unit_comment = (length_unit == "")?"": 
    "[" +length_unit + "]";
    
  Table* output = file.addTable(string (hdu), 0);
  output->addKey("LENGTHUNIT", length_unit, "Length unit for grid");
  output->addKey("MINX", getmin () [0], length_unit_comment);
  output->addKey("MINY", getmin () [1], length_unit_comment);
  output->addKey("MINZ", getmin () [2], length_unit_comment);
  output->addKey("MAXX", getmax () [0], length_unit_comment);
  output->addKey("MAXY", getmax () [1], length_unit_comment);
  output->addKey("MAXZ", getmax () [2], length_unit_comment);
  output->addKey("NX", g_.getn () [0], "");
  output->addKey("NY", g_.getn () [1], "");
  output->addKey("NZ", g_.getn () [2], "");
  
  const vec3i sub_div = get_subdiv ();
  string subdivtp ("NONE");
  bool save_subdiv = false;
  // default is no subdivision
  if (all (sub_div == vec3i (-1, -1, -1))) {
    // this means nonuniform subdivision, subdivision is saved
    // separately in another column (currently, we can't make grids like that) 
    subdivtp = "VARIABLE";
    save_subdiv = true;
  }
  else if (any (sub_div != vec3i (0, 0, 0))) {
    // uniform subdivision indicated by sub_div
    subdivtp = "UNIFORM";
  }

  output->addKey("subdivtp", subdivtp, "Type of grid subdivision");
  output->addKey("subdivx", sub_div[0], "");
  output->addKey("subdivy", sub_div[1], "");
  output->addKey("subdivz", sub_div[2], "");

  // Now go through the grid and accumulate the information in the vectors
  std::cout << "building structural information" << std::endl;
  std::vector<bool> structure;
  std::vector<vec3i> subdivision;
  g_.get_structure (structure, save_subdiv? &subdivision : 0);

  output->addColumn(Tlogical, "structure", 1,
		    "True means cell is subdivided, listed in a depth-first fashion" );
  Column& c_structure = output->column("structure");
  c_structure.write(structure, 1);

  if (save_subdiv) {
    output->addColumn(Tbyte, "subdivision", 3, "");
    Column& c_sub = output->column("subdivision");
    write (c_sub, subdivision, 1);
  }
}


/** Constructor creates grid from a grid structure FITS HDU. The
    grid_base constructor can't actually do much here.  It simply
    initializes the values of the base grid, but loading the actual
    structure has to be deferred until the cells have been created
    which happens in the constructor for the concrete grid.  That
    constructor then calls grid_base::load_recursive which creates the
    hierarchical structure. */
template<typename cell_data_type>
mcrx::grid_base<cell_data_type>::grid_base (CCfits::ExtHDU& input, 
					    const vec3d& translate_origin) 
  : parent (0)
{
  float dummy0, dummy1, dummy2; //CCfits can't read doubles
  input.readKey("minx", dummy0);
  input.readKey("miny", dummy1);
  input.readKey("minz", dummy2);
  min = vec3d (dummy0, dummy1, dummy2);
  input.readKey("maxx", dummy0);
  input.readKey("maxy", dummy1);
  input.readKey("maxz", dummy2);
  max = vec3d (dummy0, dummy1, dummy2);
  min-=translate_origin;
  max-=translate_origin;
  vec3i n;
  input.readKey("nx", n[0]);
  input.readKey("ny", n [1]);
  input.readKey("nz", n [2]);
  size = (max - min)/n;
  invsize = 1/size;
  ASSERT_ALL(size>0);

  // This is actually as much as we can do in the grid_base
  // constructor, because we need to construct the actual cells before
  // we can proceed.  The rest is performed by calling load () by the
  // concrete grid class constructor.  (Slightly weird...)
}

/** Constructor creates grid from a grid structure FITS HDU.  When
    execution enters this constructor, the grid_base constructor has
    completed. This function then looks for the refinement mechanism
    and ensures that the file is consistent with the class.  A memory
    block is allocated to hold all the subgrids, so that they are well
    localized in memory, and finally load_recursive is called to
    actually create the structure. */
template < typename cell_data_type>
mcrx:: adaptive_grid<cell_data_type>::
adaptive_grid (CCfits::ExtHDU& input, const vec3d& translate_origin) : 
  g_(input, translate_origin), subgrid_block (0)
{
  using namespace CCfits;

  string subdivtp;
  input.readKey("subdivtp",subdivtp);
  if (subdivtp == "NONE") {
    // in this case (but we have no grid like that) we're simply done
    // now
    return;
  }
  else if (subdivtp == "VARIABLE")
    // we don't know how to handle this, abort
    throw CCfits::FITS::OperationNotSupported ("Can't load variable refinement grids");
  else if (subdivtp != "UNIFORM")
    // well then what kind of a grid is it
    throw CCfits::FITS::OperationNotSupported ("Unknown grid refinement type: " +
				 subdivtp);
  else {
    // we should also make sure the divisor is appropriate for us
    vec3i subdiv;
    input.readKey("subdivx", subdiv[0]);
    input.readKey("subdivy", subdiv [1]);
    input.readKey("subdivz", subdiv [2]);
    if (any (g_.get_subdiv () != subdiv)) 
      // not correct, spew
      throw CCfits::FITS::OperationNotSupported ("Inappropriate subdivisor");
    
    // OK, now proceed to load the structure information and create it
    std::vector<bool> structure;
    Column& c_structure = input.column("structure");
    c_structure.read(structure, 1, c_structure.rows() );

    create_structure(structure);
  }
}

/** Constructor allocates cells given the grid size specified by the
    HDU keywords. */
template < typename cell_data_type,typename sub_grid_type >
mcrx:: dynamic_grid<cell_data_type,sub_grid_type>::
dynamic_grid (CCfits::ExtHDU& input, const vec3d& translate_origin) : 
  grid_base<cell_data_type> (input, translate_origin), cells ()
{
  vec3i n;
  input.readKey("nx", n[0]);
  input.readKey("ny", n [1]);
  input.readKey("nz", n [2]);
  cells.resize(n);
  
  initialize_cells ();
}

/** Recursively creates the hierarchical grid structure.  A
    depth-first traversal of the grid is done, and if a cell is
    supposed to be subdivided as indicated by the vector, the
    subdivision is done and the function calls itself in that subgrid.
    Placement new is used to make the grid structure compact in
    memory, using the supplied pointer to a preallocated memory
    block. */
template<typename cell_data_type>
void
mcrx::grid_base<cell_data_type>::
load_recursive (std::vector<bool>::const_iterator& refined,
		const std::vector<bool>::const_iterator& e,
		void*& placement )
{
  for (local_iterator c = lbegin (); c != lend (); ++c) {
    assert (refined != e);// if the file is correct this should never happen
    if (*(refined++)) {
      c->refine(placement);
      c->sub_grid()->load_recursive(refined, e, placement );
    }
  }
}

